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ABSTRACT 

We apply a recently developed scaling technique to the Millennium- XXL, one of the 
largest cosmological N-body simulations carried out to date (3 x 10^^ particles within 
a cube of volume ~ 70 Gpc^). This allows us to investigate the cosmological parameter 
dependence of the mass and evolution of haloes in the extreme high-mass tail of the 
z — Q distribution. We assume these objects to be likely hosts for the population of rare 
but ultraluminous high-redshift quasars discovered by the Sloan Digital Sky Survey. 
Haloes with a similar abundance to these quasars have a median mass of 9 x 10^^ Mq 
in the currently preferred cosmology, but do not evolve into equally extreme objects at 
z = 0. Rather, their descendants span the full range conventionally assigned to present- 
day clusters, 6 x 10^^ to 2.5 x lO^^M© for this same cosmology. The masses both at 
z — Q and at 2 = shift up or down by factors exceeding two if cosmological parameters 
are pushed to the boundaries of the range discussed in published interpretations of 
data from the WMAP satellite. The main factor determining the future growth of a 
high-mass z = 6 halo is the mean overdensity of its environment on scales of 7 to 
14 Mpc, and descendant masses can be predicted 6 to 8 times more accurately if this 
density is known than if it is not. All these features are not unique to extreme high-z 
haloes, but are generic to hierarchical growth. Finally, we find that extreme haloes 
at z = 6 typically acquired about half of their total mass in the preceding 100 Myr, 
implying very large recent accretion rates which may be related to the large black hole 
masses and high luminosities of the SDSS quasars. 

Key v^rords: cosmology:theory - large-scale structure of Universe. 



INTRODUCTION 



Bright z ^ Q quasars are extremely rare objects. Their lu- 
minosity is thought to be a result of su permassive black 
holes accreting gas at enormous rat es IjFan et al.l 120031 : 
iKollmeier et al.ll2006l : IShen et al.ll2008l ') at a time when only 
~ 10% of the mass in the Universe was in haloes where gas 
could cool efficiently and about half was still in diffuse form 
IjAngulo fc Whitel boiOal ). If these quasars are long-lived, 
then their dark matter haloes belong to an equally extreme 
tail (presumably the most massive tail) of the halo distribu- 
tion, a hypothesis supporte d by the st rong observed cluster- 
ing of high-z quasars (e.g. IShen et al. 2007). The observed 
QSO number density (about one per cubic Gigaparsec) com- 
bined with abundance matching implies that their host halo 
masses are well above 10^^ Mq at 2: = 6. These masses cor- 
respond to much more extreme peaks in the initial Gaussian 
density fluctuation field than those associated with even the 



largest galaxy clusters today. Such quasars are clearly ex- 
ceptional events in a ACDM universe, but their properties, 
together with those of the intergalactic absorption seen in 
their spectra, encrypt key inform ation about many astro- 
physical processes (e.g. iFanI 120061 ) and perhaps also about 
the background cosmological model. 

Studying the properties, environment and fate of the 
high-mass haloes in which the z ^ Q quasars may live is a 
challenging task for cosmological N-body simulations. It re- 
quires a very large computational volume to obtain a repre- 
sentative sample of extremely rare objects, sufficient spatial 
and mass resolution to resolve their structure reliably, and 
sufficient time resolution to build the merger trees needed to 
trace evolution throughout cosmic time. Previous studies in 
this area have relied on analytic mode ls calibrated using sim- 
ulations of less extreme objects (e.g. iTrenti et al.ll2008l ). or 
resim ulation techniques applied to a limited number of sys- 
tems JLi et al.ll2007l: ISiiacki et aLll2009l : iRomano-Diaz et all 
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[2OIII). In this paper, we present an extremely large N-body 
calculation, which meets the computational challenges di- 
rectly by simulating a very large volume at adequate reso- 
lution within the ACDM paradigm. 

With this simulation in hand, we study a variety of top- 
ics associated with the assembly and the future of quasar 
haloes. For example, where are today's descendants of the 
massive black holes that powered quasar s at z ~ 6? This 
question was explored by ISpringel et al.l (|2005r ) using the 
Millennium Simulation (MS). They concluded that these 
black holes would today lie at the centres of cD galaxies 
in massive galaxy clusters. However, the MS is too small to 
contain even one object like the SDSS « ~ 6 quasars, so 
this con clusion was based o n a small number of less rare 
systems. iTrenti et al.l (|2008l ) used analytic tools to extend 
these MS-based results, arguing for a much greater diver- 
sity in the present-day descendants of SDSS QSO ha l oes. A 
similar conclusion was reached by iDi Matteo et al.l (|2008l ) 
with cosmological SPH simulations. Here, we are able to use 
fully resolved N-body merger trees for haloes identified at 
the same space density as the brightest SDSS quasars. We 
confirm that these massive z = Q haloes should evolve into 
a variety of systems today. Most of their associated black 
holes should end up in the central galaxies of haloes rang- 
ing from rich groups to clusters. We show that the median 
mass of descendants approaches lO^^i\f0 but this value and 
that of early haloes themselves depend on the parameters 
assumed for the background cosmological model. 

Another issue we investigate is the large-scale environ- 
ment surrounding extreme high-redshift haloes. Such haloes 
are expected to be strongly biased towards overdense re- 
gions. Although we confirm this, we also show that there is 
considerable scatter, and many extreme haloes have envi- 
ronments of moderate overdensity. Indeed, some even have 
environments which are slightly underdense. 

Finally, we explore the assembly histories predicted for 
extreme z = Q haloes for three different assumptions about 
the parameters underlying the background ACDM cosmol- 
ogy. We measure the accretion rates onto these haloes over 
the time period immediately preceding the epoch of observa- 
tion and show that these can be extremely high. This may 
well be related to the fuelling of the extraordinarily high 
luminosities and masses measured for the 2: ~ 6 quasars. 

This paper is set out as follows. In Section 2, we present 
technical details of our simulation and of the methods we use 
to identify dark matter haloes and their evolutionary paths. 
Section 3 then presents our results, including the properties 
of the 2 = descendants of high-z quasar haloes, and their 
assembly histories at higher redshift. Our conclusions are 
presented in Section 4. 



2 NUMERICAL METHODS 

In this section we describe the numerical tools used in this 
paper. We first present our N-body simulation (Section 2.1), 
including a description of the construction of halo and sub- 
halo catalogues and of merger trees (Section 2.1.1). In Sec- 
tion 2.2 we then explain how these numerical data can be 
used to explore structure formation in cosmologies other 
than the one used to carry out the original simulation. The 



last subsection (Section 2.3) defines the samples of extreme 
haloes that we will study in the remainder of the paper. 



2.1 The MXXL N-body simulation 

Our simulation, named the Millennium-XXL or MXXL, rep- 
resents the matter content of a ACDM universe using more 
than 303 billion particles (6720^) within a comoving cube of 
side 4.1 Gpc. Our choice of cosmological parameters is iden- 
tical to t hat used in the othe r Millennium Simulations (se e 
Table 1, ISpringel et aLll2005l : iBovlan-Kolchin et al.l 120091 1 . 
implying a particle mass of 8.45 x 10^ Mq . This set of param- 
eters is inconsistent with the most recent constraints from 
observations of the Cosmic M icrowave Backg round and low- 
redshift large-scale structure JKomatsu et al. 2011 ). but, as 
discussed in Section 2.2, our results can be scaled accurately 
to any nearby cosmology, including those that are now more 
favoured. The comoving softening length of the simulation 
was e = 13.7 kpc, implying ^ 10^® effective resolution ele- 
ments in the full simulated volume. The enormous statistical 
power and dynamical range that this implies is illustrated 
in Fig. [1] which shows the 2 = density field on four differ- 
ent scales, starting from the whole box in the background, 
then zooming progressively onto the most massive halo in 
the insets. 

The initial phase-space distribution of the particles 
was set up at 2 = 63 by perturbing a glass-like distribu- 
tion (Whita 1 19961 : iBaugh et al.l 1 19951 ') using second orde r 
Lagrangian perturbation theory f2LPT IScoccimarrdll998l '). 
The use of 2LPT has several ad vantages over the more com- 
mon Zeldovich approximation (|ZerDovichlll970l ). including 
small and rapidly-decaying transients in the matter clus- 
tering and a better representation of the perturbations that 
seed extremely large dark matter haloes (|Crocce et al.ll2006l : 
iKnebe et aTll2009l : ijenkinsll2010h . The latter is particularly 
relevant for this paper. The amplitude of individual Fourier 
modes was set hierarchically as described in Jenkins (2012, 
in prep). This allows an efficient and consistent generation 
of initial conditions for the resimulation of any selected sub- 
region at, in principle, arbitrarily high mass resolution. 

Note that the MXXL is the first simulation with a large 
enough volume and a small enough particle mass to con- 
tain a representative sample of well-resolved (^ 1000 par- 
ticle) haloes which could represent the hosts of the 2 ~ 6 
SDSS quasars. Performing a simulation with these charac- 
teristics posed severe computational challenges with respect 
to raw execution time, scalability of the simulation algo- 
rithms, memory consumption and I/O performance. Thanks 
to a highly optimised simulation code and one of the largest 
supercomputers in Europe, these challenges were success- 
fully met. Specifically, the MXXL w as carried out with a 
special version of the GADGET-3 code (|Springelll2005l 'l. which 
aggressively reduced peak memory consumption at runtime, 
incorporated a number of analysis tools on the fly, and im- 
plemented a special compression of the output data. As a 
result the MXXL was completed in late summer 2010 in less 
than 3 million CPU hours (including postprocessing calcula- 
tions), using 30 Tb of RAM, and 12,228 cores of the Juropa 
cluster at the Jiilich Superc omputer Center in Germany. We 
refer to I Angulo et al.l (|2012l ) for more details about the sim- 
ulation. 
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Figure 1. Images of the matter density field in the Millennium-XXL focusing on the most massive halo present in the simulation at 
2 = 0. Each inset zooms by a factor of 8 from the previous one; the side-length varies from 4.1 Gpc down to 8.1 Mpc. The intensity 
of each pixel is proportional to the logarithm of the dark matter density projected through a 25 Mpc thick slab. This simulation has a 
dynamic range of 3 X 10^ on each spatial dimension, simultaneously resolving the internal structure of collapsed objects and the large-scale 
quasi-linear fluctuations in a ACDM universe. 



2.1.1 Halo and Subhalo catalogues 

We identified self-bound halo/subhalo structures through- 
out the MXXL at the same 64 redshifts used in the MS and 
MS-II. This output frequency (roughly equally spaced in 
time by 300 Myr for z <2, and by 100 Myr at z ~ 6) allows 
us to build detailed merger trees. At each out put time, we 
first apply a Friends-of-Friends (FoF) algorithm (jDavis et al.l 
1 1985 ), with a linking length of 0.2 times the mean interpar- 
ticle separation to build a FoF group catalogue down to a 
limit of 20 particles. We then use a m emory-efficient imple - 
mentation of the SUBFIND algorithm (|Springel et al.l 120011 ) 
to identify self-bound substructures within each FoF group 
down to a limit of 15 particles. These calculations were per- 



formed on-the-fly during the N-body calculation, so that it 
was not necessary to store the particle data at all output 
times. This significantly reduced the I/O and storage re- 
quirements of the simulation. 

Summing over all output times, there are 2.5 x 10^" FoF 
groups in the MXXL with more than 20 particles. At 2 = 6 
there are 3.7 x 10^ such groups and 6.5 x 10* at z; = 0. The 
most massive FoF group at 2: ~ 6 contains 3,285 particles 
and 4 substructures. This is about 300 times less massive 
than the biggest halo at the z — snapshot, which contains 
1,062,232 particles and 688 substructures with more than 15 
particles. 

Finally we built "merger trees" similar to those de- 
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mdn 



Box 



f^n 



Qh 



0"8 



MXXL(Ml) 0.95 
M3 1.64 

M7 1.26 



4167 
5091 
4467 



0.250 
0.238 
0.272 



0.045 0.9 
0.0416 0.761 
0.0416 0.807 



Table 1. Parameters of the original MXXL simulation and of the 
scaled versions used to represent other cosmologies. The columns 
are as follows: (1) the name of the simulation; (2) the mass of 
a dark matter particle in units of lO^'' Mq ; (3) the side of the 
computational box in units of Mpc; (4) the total matter density; 
(5) the baryon density; (6) The linear fluctuation amplitude at 
2 = 0. In all cases the primordial spectral index is Us = 1, the 
Hubble constant at 2 = is Hq = 73 km/s/Mpc, and the dark 
energy is assumed to be a cosmological constant. 



scribed in lSpringel et al.l ([20051) in order to follow the evolu- 
tion of halo/subhalo structure in detail. For every subhalo in 
our catalogues we define a pointer to a unique descendant in 
the subsequent snapshot by locating the subhalo containing 
the greatest number of its 15 most bound particle^ These 
pointers are used to create a tree-like data structure, which 
represents the full assembly history, the current substruc- 
ture, and the future evolution of every halo. In particular 
this allows us to map out the formation histories of our pu- 
tative z = Q quasar hosts, and to follow their later evolution 
down to 2 = 0. 



2.2 Exploring structure formation in other 
cosmologies 

The MXXL was an extremely expensive numerical simula- 
tion and so could only be carried out once for a specific set of 
cosmological p arameters. However, the rescaling method of 
lAngulo fc White (2010h ) allows us to use the simulation to 
analyse any neighbouring cosmology with Gaussian initial 
fluctuations, and in this paper we will show results not only 
for the original MS cosmology but also for two other ver- 
sions of the standard ACDM cosmology. The accuracy of the 
rescaling scheme is remarkably high if it is applied carefully; 
masses of individual objects are reproduc ed to better than 
10% and positions to be tter than 100 kpc (jAngulo fc White! 
uOlObI : iRuiz et al.ll2011h . In the following, we briefly recap 
the main features of the method. 

First, consider a "target" cosmological model a,t z = zb 
which we seek to match using the results of an "origi- 
nal" cosmological simulation of side-length La (in units of 
h~^ Mpc). The heart of the method is to find a length trans- 
formation. La — >■ Lb ~ s La, and a relabelling of the time 
variable za ^ zb, by requiring that the variance of the lin- 
ear density field in the target cosmology, a%{R,ZB), over 
the range [Ri, R2] is as close as possible to that of the origi- 
nal cosmology, (T^(i?, z;a), over the range [s~^Ri,s~^R2] at 
redshift za- Thus, we minimise: 



R2 



Ri 



dR 



'oI{R)Db(zb) ~ g\{s-^R)Da{za) 



(1) 



If two subhaloes contain the same number of these particles, 
we choose the one with the largest total binding energy. 



over a and za, where D{z) is the linear growth factor in 
units of its present-day value. 

In the Press-Schechter theory ()Press fc Schechtejll974h 
the halo mass function is determined by the linear vari- 
ance of the underlying dark matter field, thus Eq. ([T| also 
minimises the difference between the halo mass functions 
in the target and original cosmologies over the mass range 
[M{R2),M{Ri)]. We usually take M{R2) to be the mass of 
the largest halo in the simulation at the lowest redshift of 
interest {z = here) and M{Ri) to be that of the least 
massive resolved halo. 

As a result, the original box size will be expanded by 
a factor s, and the output at redshift za will represent red- 
shift zb in the target cosmology. Redshifts in the target 
cosmology (z) corresponding to the redshifts (z) of stored 
data in the original simulation are then determined implic- 
itly by Db{z') = [Da{z)/Da{za)]Db{zb). The mass of a 
simulation particle (in units of Mq) in the target cosmology 
is m_g = {Q.m,BH%/Q.m,AHA) s^ ruA, where Q.m,x and Hx 
[X = A ox B) are the dimensionless total matter densities 
and Hubble parameters of the two cosmol ogies. 

Th e second part of the algorithm of lAngulo fc White! 
(|2010b|) corrects differences in power spectrum shape on 
large scales between the original and target cosmologies 
by altering the amplitude of quasi-linear modes using the 
Zel'dovich approximation. In this paper we do not make this 
correction since we are interested in the internal structure, 
abundance and evolution of massive haloes rather than in 
their spatial distribution. 

With this technique we have created two additional halo 
catalogues in alternative cosmologies which we denote M3 
and M7. T hese have cosmologi cal paramete rs motivated by 
the 3-year (jSpergel et al.!!2007! ) and 7- year (JKomatsu et al.! 
I^OII) analyses of data from the WMAP satellite. The main 
feature of these models are values for ag, which are lower 
than in the MXXL and the other Millennium Simulations, 
and different values for Q,m (see Table 1). The corresponding 
length scalings are s = 1.222 and s = 1.072 for MS and M7, 
respectively. The z = 0.623 and z — 0.319 outputs of the 
MXXL represent z' = in the M3 and M7 cosmologies. 



2.3 QSO haloes 

In this paper we will assume that the QSO luminosity is 
a monotonically-increasing function of the FoF host halo 
mass at any given time, and that there is a duty cycle that 
is independent of the halo mass. Models with these charac- 
teristics appear to be preferred by clustering analyses (e.g 
I White et aLll2008l : lBo"noh et al.ll2010! '). but we note that they 
are not the only possibility. In physically motivated models 
the QSO host depends on the details of the gal axy formation 
model and in particula r on AGN feedback (JMarulli et al.! 
I2OO8! : JBonoh et al.1l2009l : JFanidakis et"alll201l! . l20li v 

Therefore, we consider halo samples limited by FoF 
mass at two different thresholds corresponding to two differ- 
ent abundances at z ^ 6 which we keep constant across our 
three cosmologies, "Long-lived QSO haloes" (LLQ haloes) 
have a comoving number density of n = 0.4 Gpc~^ , requiring 
FoF masses above [15.4, 9, 5] x 10^^ A'/q for the Ml, M7 and 
M3 cosmologies, respectively. The number density of this 
sample (which is well below the limit that could be probed 
by the MS) matches that of the extremely bright quasars 
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observed a t z > 5 in the SDSS (n = 0.6 Gpc"^, iFan et all 
12001 I2OO6I '). Thus, they mimic the assembly history of SDSS 
QSOs if they have a 100% duty cycle, i.e. if they shine con- 
stantly at their full brightness. 

Our second sample, "Short lived QSO haloes" (SLQ 
haloes) is selected at 30 times higher abundance, i.e. n = 
11.6 Gpc~^. This corresponds to minimum FoF halo masses 
of [7.0, 4.2, 2.0] X 10^^ Mq for the Ml, M7 and M3 cosmolo- 
gies. This sample could be regarded as representing the hosts 
of the high-redshift SDSS quasars if these have a 1/30 duty 
cycle, i.e. each object shines at full brightness only 1/30 of 
the time. 

Due to our length scaling, the total number of haloes 
in our two samples varies by factors of s^ between the three 
cosmologies; the SLQ samples contain 810, 997 and 1478 
objects for the Ml, M7 and M3 cases, respectively (note 
that we expect four objects above this mass threshold in 
the MS which, in fact, contained only two), whereas the 
LLQ samples contain 27, 34, and 50 haloes. The redshift 
of the samples also differs between cosmologies because the 
MXXL data were stored only at a discrete set of times: we 
use 2 = 6.18 (Ml), 2 = 6.19 (M3) and z = 5.96 (M7) for the 
three cases. 

Finally, we note that alternative sample definitions, for 
example, using virial masses or circular velocities rather 
than FoF masses, do not produce significant differences in 
our results. 70% of the 1000 most massive FoF haloes rank 
within the 1800 most massive haloes according to M200, and 
within the 4000 according to peak circular velocity. 



3 RESULTS AND DISCUSSION 

Using the numerical tools and halo samples presented in the 
previous section, we now study the evolution of the host 
haloes of z ~ 6 quasars. We look into the type of objects 
they turn into at the present day (Section 3.1) and how the 
masses of these descendants are related to the environments 
of their 2 ~ 6 progenitors (Section 3.3). In addition, we ex- 
plore the growth rates of our quasar host haloes both prior 
(Section 3.4) and subsequent (Section 3.2) to their identifi- 
cation at 2: ~ 6. 



3.1 The mass and fate of 2 ~ 6 quasar hosts 

Since we are assuming that high-redshift quasars live in the 
most massive objects present at that time, it seems plau- 
sible that their descendants will lie at the centres of the 
most massive haloes at any later time, in particular, in 
the central ga l axies of large galaxy clus ters today (see e.g. 
ISpringel etall l2005l : ICapak et all bOllh . However, in this 
section we show that this is not always true. The descen- 
dants of QSOs can be found in haloes spanning a factor of 
30 in mass, and in a few cases, they are not even located in 
the central object of this halo, but in an orbiting subhalo. 

The upper panel of Fig. [2] shows differential mass func- 
tions for the two 2 ~ 6 QSO host halo samples (defined in 
Section 12. 3p and for the three cosmological models of Ta- 
ble 1, which differ primarily in their value of ag, and Q.m- By 
construction, these samples consist of very massive haloes. 
Their mean mass ranges from 2.6 to 8 x 10^^ Mq for the 
SLQ samples, and from 5.2 to 15.2 x lO'^^ Mq for the sparser 
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Figure 2. Comparison of the FoF mass distributions of our 
SLQ host halo samples at 2 ~ 6 to those of their descendants 
at z = 0. The top panel shows the differential number density of 
haloes for the three cosmologies of Table 1, corresponding to the 
parameters preferred by analyses of one, three and seven years of 
data from the WMAP satellite. On each histogram a thick vertical 
line indicates the mass threshold for the corresponding LLQ host 
halo sample. The bottom panel shows the distributions of FoF 
masses of the z = d escen dants of these h aloes. Predictions of the 
I Jenkins etaP l|200lh and lAngulo et al] ll2012l) fitting formulae , 
and of merger trees constructed using the lParkinson et al.l ||2008J) 
algorithm are shown in the top and bottom panels, respectively. 
Vertical arrows in the bottom panel indicate the median masses 
of the descendants of SLQ haloes (n = 30Gpc~ ) and of LLQ 
haloes (n = 1 Gpc~ ). Note that haloes of similar mass at z ~ 6 
end up in haloes with a wide range of masses at z = 0. 



LLQ samples - for both samples these mean masses increase 
smoothly with erg, as expected. Within each sample more 
than than 99% of all haloes lie within a factor of two in 
mass of the threshold for inclusion. This is a consequence 
of the exponentially falling high-mass tail of the halo mass 
function. 

The magnitude of the shifts between the three cos- 
mologies are well de scribed by the I Jenkins et al.l (|200ll ') and 
lAngulo et al.l ([201^) fitting formulae, which we display as 
solid or dashed lines in Fig. [2] However, at this redshift these 
formulae overpredict the number of haloes of a given mass by 
a factor of two to three. This reflects the fact that these for- 
mulae were calibrated using simulations and redshifts where 
the most massive haloes were much less extreme than those 
considered here, but the disagreement might also be a con- 
sequence of the non-universal beha viour of the halo mass 
function (see e.g [Tinker et al.ll2008l ). 
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We identify the z = descendant of each halo in our 
samples as the object that contains the majority of its 15 
most bound particles. Although in principle it is possible 
that a different structure contains most of the mass of our 
high-2 haloes, following the innermost particles should rep- 
resent well the fate of a hypothetical black hole sitting at 
the centre of the halo. The bottom panel of Fig. [2] displays 
the FoF mass distributions of these descendant haloes. For 
comparison, we also present p redictions made using t he ana- 
lytic merger tree algorithm of I Parkinson et al.l l|2008i ) which 
is based on the Extended Press-Schechter (EPS) formalism 
and calibrated using merger trees extracted from the Mil- 
lennium Simulation. 

The median FoF mass of the descendants of our SLQ 
haloes varies by a factor of just 1.7 across our three cos- 
mologies (see the lower set of coloured arrows in Fig. [2]). 
This is significantly smaller than the spread of a factor of 
3.1 in the the initial median masses. The typical descendant 
mass is M ~ 5.6 x 10^"* Mq, corresponding to a moderately 
rich 2 = cluster. The median masses of the descendants 
of the sparser LLQ haloes are more massive by a factor of 
about 1.5 and also vary slightly more with ag (see the upper 
set of coloured arrows in Fig. [Sjl . The median initial masses 
of the LLQ haloes are typically a factor of 1.9 larger than 
those of the SLQ haloes, so both within a single cosmology 
and between cosmologies the evolution is convergent in the 
sense that descendant masses are more similar than those of 
the original 2; ~ 6 objects. This is probably a result of QSO 
haloes, in all cases, descending into much less extreme peaks, 
where the differences among cosmologies are smaller. For ex- 
ample, the mass function for the Ml and M7 cosmologies are 
almost identical at z = for masses below M ~ 10^* Mq . 

In all cases, there is a large spread in the masses of the 
descendants, as shown by the bottom panel of Fig. [S] This 
distribution at z = can be well described by a log normal 
with o — 0.36 dex. The most massive tail indeed corresponds 
to rare, high-mass clusters, but the least massive one cor- 
responds to galaxy groups. The scatter is slightly smaller 
for the LLQ sample. We also would like to note that every 
descendant of an LLQ halo is the dominant structure of its 
z = group, but in the Ml cosmology 39 out of the 810 
descendants of SLQ haloes are satellite subhaloes. For the 
M3 and M7 cases the corresponding numbers are similar. As 
we will see in the next subsection, all this reflects the variety 
of formation histories among dark matter haloes of similar 
mass. 

Our distributions are in goo d qualitative agree ment 
with the EPS-based calculations of iTrenti et~all l|2008h . Us- 
ing the same Ml cosmology, these authors found 68% of 
the descendants of a sample analogous to our LLQ to have 
masses in the range 2.5 to 12.2 x lO'^* Mq with a median 
of 5.6 X 10^'^ Mq. For this cosmology, our results show a 
similar scatter, but with a median which is 35% larger. As 
can be seen in Fig. [2] similar s catter is fo und if we use the 
merger tree algorithm of Parkinson et al.l ( 2008 ) to predict 
the masses of the descendants of our SLQ haloes, but this al- 
gorithm also predicts a median mass which is too small by a 
factor of about 1.5. This suggests that the results from EPS 
should not be extrapolated to objects much more extreme 
than those for which they were calibrated. 
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Figure 3. Mass growth of high-redshift quasar host haloes in 
the Ml cosmology (i.e. in the original MXXL). Blue regions show 
mass growth for the SLQ sample, whereas brown regions show 
growth for the sparser and more massive LLQ sample. In each 
case, the dark and light regions enclose 68% and 86% of the tra- 
jectories. Individual lines highlight particular trajectories. The 
thick solid line is the growth history that the largest halo at z = 6 
would have if its mass was equal to that of the most massive halo 
in the simulation at all later redshifts. In contrast, the dotted line 
shows the actual growth history of this halo. Finally, the dashed 
line displays the growth of the lowest mass halo in the SLQ sam- 
ple which, by chance, is one of the fastest growing of all haloes. 



3.2 The journey to 2 = 

We now look in more detail at the paths connecting our 
samples of z '^ 6 quasar host haloes to their present-day de- 
scendants, concentrating on results in the unsealed MXXL, 
i.e. in the original Ml cosmology. In Fig. |3] we plot mass 
growth, defined as the ratio of descendant mass at each 
redshift to initial mass at 2 ~ 6. Light blue and brown 
regions indicate 86% of the trajectories for SLQ and LLQ 
samples, respectively. Darker regions of each colour outline 
the regions containing 68% of the trajectories. Growth his- 
tories seem to be remarkably similar in the two samples and 
to be faster than exponential, described approximately by 

0.21 X (6.19 - 2)^ ^ This seems 



logio A'/(2)/M(z = 6.19) 
independent of initial mass at 2 = 6, at least over the rela- 
tively restricted range of (high) masses considered here. Of 
course, this formula only describes the typical behaviour, 
and very different growth histories occur for haloes of simi- 
lar initial mass. Some massive 2 = 6 haloes have grown orfly 
by a factor of 20 to 30 by 2 = - much less than the change 
of the characteristic nonlinear mass-scale M« over the same 
redshift range - while others have increased their mass by 
factors of 200 to 300. 

The spread in these trajectories increases substantially 
with decreasing redshift. At 2 = 3 the rms scatter in the log 
of the fractional growth is 0.168, at 2 = 1 it is 0.271 and 
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Figure 4. Fraction of descendants of SLQ haloes that are among 
the 11 (soUd), 112 (dotted) or 1120 (dashed) most massive haloes 
per Gpc at each later time (top panel), or that have a mass 
above 2.8 X lO" Mq, 1.4 X IQl-* Mq or 1.4 X lO^^ Mq (bottom 
panel). In the top panel, green, red and black lines correspond to 
the results for the Ml, M7 and M3 cosmologies. 



at 2 = it is 0.33, slightly larger than the initial separa- 
tion in median mass between our SLQ and LLQ samples. 
Although we do not display them, the mass accretion histo- 
ries of haloes in other cosmologies are very similar. This is, 
of course, expected. 

Fig. [3] also shows the mass growth for the most massive 
(dotted) and least massive (dashed) haloes in our SLQ sam- 
ple. By chance, the least massive halo is one of the fastest 
growing, with a fractional growth rate faster than that of the 
most massive halo at all redshifts. Indeed, this "low-mass" 
halo grows into a 2 x 10^^ Mq object by z = , whereas the 
present-day descendant of the initially most massive halo is 
actually slightly smaller (1.8 x IQ^^ Mq). Neither of these 
haloes is in the extreme tail at 2 = 0, ranking 2,224-th and 
3,459-th in mass among MXXL haloes, and being four times 
less massive than the most extreme object. Curiously, none 
of the 20 most massive z = Q haloes in the MXXL has a 
progenitor in the SLQ sample. 

Finally, the thick solid line in Fig. |3] indicates the mass 
growth that the most massive halo at 2 = 6 would have to 
have in order to be the most massive halo at all redshifts. 
This is close to the actual trajectory of the largest halo until 
2 ~ 3, but at later times, other objects take over the top 
spots. 

All these examples are not exclusive to extremely mas- 
sive haloes at high redshifts, but are a generic illustration 
that the most massive haloes at any given epoch will no 
longer be the most massive haloes at a later time. This is a 
consequence of the diverse assembly histories of haloes in a 
hierarchical universe. 

We explore this behaviour directly in Fig. |4l The top 
panel indicates the fraction of SLQ halo descendants whose 



Figure 5. Correlation between the overdensity in which a QSO 
halo sits at 2 = 6 and the mass of its descendant at 2 = 0. The 
horizontal axis shows Gaussian-smoothed overdensity in units of 
its rms value (T2.5 = 0.193, 0-5 = 0.125, cio = 0.071 and cr20 = 
0.036, where the suffix indicates the 1-D smoothing radius of the 
Gaussian in units of Mpc comoving. Different colours refer to 
different smoothing scales as indicated in the figure. Small and 
large points show the overdensities of the environments of SLQ 
and LLQ haloes, respectively. 



masses place them above three abundance thresholds. The 
bottom panel shows a complementary picture, showing the 
fraction of SLQ haloes that rise above various mass thresh- 
olds at later times. 

At redshift 3, all descendants of SLQ haloes have masses 
above 2 x lO'^^ Mq and rank among the 3000 most massive 
haloes per cubic Gigaparsec. (Recall that initially these ob- 
jects are defined as the high-mass tail of the halo distribu- 
tion with an abundance of 30 per cubic Gpc.) In contrast, 
only 5% are in haloes with M > 1.4 x 10^* Mq and only 
20% still rank among the 30 most massive haloes per cubic 
Gigaparsec. With time, the SLQ descendants fall further 
and further behind the most massive haloes in the MXXL. 
By the present day, only 2 — 5% (depending on cosmol- 
ogy) would still be included in a mass-limited sample with 
n — 11.6 Gpc~^ and about 50 — 70% would be included in a 
sample with 100 times greater abundance. 



3.3 The environment of QSO haloes 

Why do some high-redshift haloes keep growing rapidly un- 
til the present, whereas others appear to shut down their 
accretion? Is this a random process or can it be related to 
some halo property at 2 ~ 6? 

FigOshows how the masses of the 2 = descendants of 
our SLQ and LLQ haloes correlate with their environment 
densities. The overdensity surrounding each high-redshift 
halo is computed on a variety of scales by mapping the un- 
derlying dark matter distribution onto a 2048'^ grid and then 
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convolving this density field with a Gaussian filter of differ- 
ent sizes. 

Clearly, most of the QSO haloes live in overdense re- 
gions, but there is considerable diversity in how extreme 
these regions are. On small scales, all live in at least a 3(t 
region, but la is the typical overdensity of an SLQ halo and 
10(7 that of an LLQ halo. With increasing smoothing, our 
quasar host haloes are found in progressively less extreme 
regions and the separation between the SLQ and LLQ haloes 
decreases. For a smoothing of 28 Mpc the typical SLQ halo 
lives in a region of overdensity 1.4cr, and 10% are located in 
regions of below-average density. Thus, our results suggest 
that it might be possible to find a quasar in the middle of 
a 14 — 28 Mpc underdense region, even if quasars reside in 
the most massive haloes at 2: ~ 6. 

Fig. [5] also shows that there is clearly a strong correla- 
tion between the overdensity around a high-2 halo and the 
mass of its z = descendant. In fact, this correlation is 
stronger than between the mass of the z = Q progenitor and 
that of its 2; = descendant (or with any other property 
in our catalogues). For example, if we consider the environ- 
ment at 14 Mpc, haloes that live in < la regions end up in 
haloes of M ~ 2.8 x 10^"' Mq, whereas those found in > 4cr 
regions typically end up in ten times more massive haloes 
(2.8 X 10^^ Mq). The scatter in descendant mass at a fixed 
overdensity is aiogM = [0.061, 0.044, 0.045, 0.063] for the 
fields smoothed on [3.5, 7, 14, 28] Mpc scales, respectively. 
These figures are to be compared with a scatter of 0.36 dex 
for the sample as a whole. Thus, if the environment density 
surrounding a quasar is known, then the mass of its z = Q 
descendant can be predicted 6 — 8 times more accurately 
than if it is not. 

These results are easy to understand, since the masses 
of the z = Q descendants correspond to the mass contained 
within a sphere of radius about 7—14 Mpc, and an object 
can only form by « = if its material is already overdense 
by a factor of about 1.68/D+(2 = 60= 0.32 at 2 = 6. 
Our findings also warn against a naive connection between 
objects at different redshifts - the linkage depends not only 
on the actual properties of the objects, but also on their 
environment. 



3.4 Prior accretion histories 

Our 2 ~ 6 quasar host haloes are the most massive objects 
present at that time and so, by definition, are the objects 
which had the highest mean mass growth rates averaged over 
previous epochs. It is these extreme growth rates which must 
supply the baryons needed to build up the supermassive 
black holes and to fuel the extraordinarily luminous quasars 
which assume to lie in their cores. In this final section, we 
examine when and how fast our LLQ and SLQ haloes achieve 
their extreme masses. 

We define an accretion history for each halo in our sam- 
ples from the growth in mass along the main branch of its 
assembly tree. This branch is defined by stepping back in 
time from the 2^6, object, selecting the most massive 
progenitor of the current main branch object as the main 
branch object at the immediately preceding step. Note that 



Z)-l-(z) is the growth factor in units of its present-day value. 
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Figure 6. The distribution of formation times for haloes in our 
SLQ sample. We define the formation time of a halo as the epoch 
when the main progenitor of a halo had 50% or 25% of the final 
halo mass, which we display in the top and bottom panels respec- 
tively. These values are given in terms of the age of the Universe 
at the time the SQL sample is identified and correspond to 0.92, 
0.93 and 0.94 Gyr for the Ml, M3 and M7 cosmologies, respec- 
tively. Dotted lines in both panels indicate the median values for 
each sample. 



with this algorithm the main progenitor at 2: = 10, say, is 
not necessarily the most massive of all the z = 10 progen- 
itors. In fact, only about ~ 40% of our SLQ haloes have 
an identifiable main progenitor at 2 '^ 10 (i.e. with a mass 
above the resolution limit, 20 particles or 1.2 x 10^^ Mq) 
even though the MXXL contains 125,000 identifiable haloes 
at this time. Conversely, of the 100 most massive MXXL 
haloes ai z — 10, only 30 have a descendant among our SLQ 
sample at 2 ^ 6. This is a further manifestation of the broad 
scatter in assembly histories found in hierarchical clustering 
models. It su ggests that any lumin ous objects discovered at 



10 (e.g. iBouwens et al.l l2011f ) are unlikely to be the 



progenitor of the 2 ~ 6 SDSS population, just as these are 
not typically the progenitors of the most massive cluster cD 
galaxies today. 

A consequence of these statistics is that SLQ haloes 
typically accrete most of their mass in a relatively short pe- 
riod before they are identified. This is illustrated explicitly 
in Fig. [S] which shows histograms of the 25 and 50% growth 
times of haloes in each of our three SLQ samples. These 
are defined for each halo as the times since it had a quarter 
and a half of its final mass, and they are given in units of 
the age of the universe, in, at the time the samples were 
defined. Median values are ipx/iH = 0.89 and 0.78, respec- 
tively for our two definitions of formation time, and they are 
almost independent of the cosmological model. These values 
correspond to roughly 100 and 200 Myr prior 2 ~ 6. 

Recent accretion rates are clearly substantially larger 
than the mean value required to grow the halo in the Hub- 
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ble time. Median accretion rates for the last half of halo 
growth are [9.7, 8.9, 7.6] x 10^ Mq yr"^ for Ml, M7 and M3 
samples, respectively. The growth times of our LLQ haloes 
are similarly distributed to those shown in Fig. |6] and, as a 
result, their median accretion rates are two to three times 
higher. In all three cosmologies, the growth rates we find 
appear to be comfortably large enough to fuel even quasars 
as bright as the SDSS objects at z --^ 6, provided, of course 
that the associated baryons are able to shed most of their 
angular momentum and reach the central regions despite the 
tremendous luminosity being generated there. 



4 CONCLUSIONS 

In this paper we have combined the largest high-resolution 
cosmological simulation to date with a scaling technique 
which allows a simulation to represent structure growth in 
cosmologies other than that in which it was originally car- 
ried out. This allows us to explore the properties and the 
evolution of extremely massive haloes that might host 2: ~ 6 
quasars in three cosmologies with parameters spanning the 
observationally allowed range. 

We found significant differences in the growth of such 
haloes subsequent to their identification at z ~ 6. Some in- 
crease their mass by a factor of 200 by 2 = whereas others 
grow only by a factor of 10. As a result, the descendants 
of bright high-redshift quasars are inferred to live in haloes 
with a wide range of halo masses today. The median descen- 
dant mass of haloes in a mass-limited sample with space 
density 11.6 Gpc^ at z ~ 6 is 5.7 x lO" Mq for a WMAP7- 
like cosmology, while more massive objects with 30 times 
lower abundance, thus matching the directly observed num- 
ber density of luminous SDSS quasars, end up in haloes with 
median mass about a factor of two higher. In both samples 
descendants spread in mass by a factor of several above and 
below this median. Conversely, in this same cosmology, only 
4% of present-day haloes with mass above 2.8 x 10^^ Af©, 
corresponding to space density 11.6 Gpc^, have a progeni- 
tor at 2 ~ 6 with mass above 7.1 x 10^^ Mq and so would 
be considered a potential quasar host at this same abun- 
dance. These figures change only slightly for the other two 
cosmologies we consider. 

Another aspect of the same effect is that haloes ranked 
among the most massive at a given time, will gradually oc- 
cupy lower positions and other haloes, initially less mas- 
sive, will eventually take over the top positions. The dissim- 
ilar mass growth is also expected to influence the galaxies 
that would form in these haloes, which violates the core as- 
sumption of HOD modelling and simple abundance match- 
ing techniques. In addition, we find that the best way to 
predict the later growth of z = 6 haloes is to look at their 
local environment on 14 Mpc scales, which correlates with 
the halo mass a,t z = Q much more strongly than the actual 
mass at 2 ~ 6. We emphasise that the behaviour we describe 
in the paper is not restricted to z ~ 6 haloes, but it is a gen- 
eral feature expected in hierarchical growth, where the ini- 
tial amplitudes of different Fourier modes are independent of 
each other. This behaviour is a example of a general property 
of certain mathematical distributions known as "regression 
to the mean" , which describes the migration of an extreme 



sample to a less extreme one at a later timelj The "regres- 
sion to the mean" can only be avoided if there were a perfect 
and monotonically increasing relation between the mass of a 
halo at that of its descendant. In this case the rank order of 
haloes by mass is perfectly preserved. Thus, extreme haloes 
at, for instance, z = Q beget equally extreme haloes at 2 = 0. 
However, we have shown that this situation is unlikely for 
hierarchical structure formation. 

Finally, we explored the assembly of QSO haloes prior 
to their identification at z ~ 6, finding one of the fastest 
accretion rates ever seen in simulated objects: a median of 
about 8.6 X 10'^ Mq yr~^ (but up to a factor 2 larger) for the 
100 Myrs preceding identification at 2 ~ 6, almost indepen- 
dent of cosmological model. This accreted mass appears to 
provide sufficient gas to fuel the bright quasars observed in 
the SDSS. 
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